Recoiling from a kick in the head-on collision of spinning black holes 
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Recoil "kicks" induced by gravitational radiation are expected in the inspiral and merger of 
black holes. Recently the numerical relativity community has begun to measure the significant 
kicks found when both unequal masses and spins are considered. Because understanding the cause 
and magnitude of each component of this kick may be complicated in inspiral simulations, we 
consider these effects in the context of a simple test problem. We study recoils from collisions of 
binaries with initially head-on trajectories, starting with the simplest case of equal masses with no 
spin and then adding spin and varying the mass ratio, both separately and jointly. We find spin- 
induced recoils to be significant relative to unequal-mass recoils even in head-on configurations. 
Additionally, it appears that the scaling of transverse kicks with spins is consistent with post- 
Newtonian theory, even though the kick is generated in the nonlinear merger interaction, where 
post-Newtonian theory should not apply. This suggests that a simple heuristic description might be 
effective in the estimation of spin-kicks. 

PACS numbers: 04.25.Dm, 04.25.Nx, 04.30.Db, 04.70.Bw, 95.30.Sf, 97.60.Lf 



I. INTRODUCTION 

The coalescence of spinning black holes in a binary 
system is expected to occur throughout the universe, 
on scales ranging from stellar black holes formed as the 
end-products of stellar evolution to supermassive black 
holes that lurk at the centers of galaxies. The final 
merger of such systems will produce an intense burst 
of gravitational radiation; if this radiation is emitted 
asymmetrically, as in the case of unequal masses and 
spins, the resulting remnant black hole will experience 
a recoil kick. The magnitude of this kick is impor- 
tant in a variety of astrophysical situations, such as 
the cosmological evolution of supermassive black holes 
[3, S IE 11 IE IE 0, II and the growth and retention of 
intermediate-mass black holes in dense stellar clusters 
[E EE HI El EE El El, and it also affects the ex- 
pected rates of black hole mergers for gravitational wave 
detectors [IE]. Given the importance of recoil kicks in 
astrophysics, there have been numerous analytic studies 
of this phenomenon [H El EE HE EH, 113, Hi, EH- 
However, since nearly all of the recoil occurs in the regime 
of strong gravitational fields, numerical relativity simula- 
tions are essential to obtain accurate calculations of the 
kick. 

Recent breakthroughs in binary black hole simulations 

he si mm 

now allow extensive studies of equal mass 
nonspinning black hole mergers HE EH, EH . Calculations 
of the kick resulting from mergers of nonspinning black 
holes with unequal masses have been carried out for mass 
ratios q = mi/mg in the range q = 1 to q = 0.25 HE 
EH EH i with the maximum kick velocity estimated to be 
~ 175km/s for mass ratio q ~ 0.36 [35|. In addition, 
several simulations of mergers of equal mass, spinning 



black holes have been carried out [IE El El] . 

Most recently, several papers have appeared that ad- 
dress the kicks obtained from the inspiral of spinning 
binaries of equal HE E3] and unequal masses [4lj| . These 
papers find considerable kicks resulting from the addi- 
tion of symmetry-breaking spins perpendicular to the or- 
bital plane, leading to a total kick speed of ~ 440km/s. 
Moving to spins initially in the orbital plane, [13] sug- 
gested a "supcrkick" orbital configuration; this was in- 
vestigated by Gonzalez et al. [43[, who found an asso- 
ciated kick of ~ 2500km/s. Extrapolated to maximally 
spinning pre-merger holes, this may lead to a kick as high 
as - 4000km/s ga]. 

Nevertheless the complexity of the binary orbital evo- 
lution may obscure details such as the direction of the 
final kick, and its dependence on mass ratio and spin. 
With this in mind, we present here a study of a sim- 
pler problem, which may be seen as an approximation to 
the final plunge to merger. We chose the head-on case 
as a model problem to isolate kick effects from the or- 
bital inspiral motion. Although lacking in astrophysical 
likelihood, our head-on investigations have the advantage 
of removing directional ambiguity in the kicks produced. 
In particular, we can readily test the leading-order post- 
Newtonian (PN) prediction that the spin and mass-ratio 
contributions to the kick should be orthogonal. 

In our investigation of the recoil kicks produced by the 
merger of spinning black holes in head-on collisions, we 
vary both the black hole spins and mass ratio. We find 
definite transverse kicks from the merger of equal-mass 
holes with spin. Furthermore, the total kick momentum 
imparted appears to scale roughly with the sum of the 
black hole spin parameters, a± + a?. For the cases inves- 
tigated, these kicks yield transverse kick velocities in the 
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range ~ 15 — 30km/s. We have also seen longitudinal 
(along the line of motion of the holes) kicks, due to a 
non-unity mass ratio. These kicks are much more mod- 
est in magnitude, with velocities of ~ (2 — 5)km/s for the 
mild mass ratio q = 2/3 chosen. 

This work is carried out using the moving puncture 
method [H [H, 0, [H, [H, E3 • We describe our initial 
data in Sec. Inland our numerical methodology in Sec. lIIII 
Code calibration and testing is presented in Sec. IIV1 Wc 
present our results in Sec. [Viand conclude with a discus- 
sion in Sec. [VlJ 

II. INITIAL DATA 

We set up initial data for binary black holes repre- 
sented as "punctures" [4^ . The metric on the initial 
spacelike slice is written in the form gij = ip^Sij, where 
i,j = 1,2,3, with conformal factor ip = -0BL + u. The 
singular part of the conformal factor takes the form 
V-'bl = 1 + X1a=i TO A,p/2|r — FaI, where the A th black 
hole has mass parameter (or "puncture mass") mA, P and 
is located at coordinate position r*A- The nonsingular 
function u is calculated by solving the Hamiltonian con- 
straint equation using the second-order-convcrgcnt ellip- 
tic solver AMRMG 0]. 

The black holes start out at rest. We use the Bowcn- 
York (BY) [50] form of the extrinsic curvature to incor- 
porate black hole spin: 
3 

K A,ij = V^ 2 — [e ktm S^n A n Aj + e kjm SJn A nAi] , (1) 

where ua is the unit vector in the direction of increasing 
rA and Sa is the spin angular momentum of the A th 
black hole. In all cases we take 5a to be aligned with the 
z axis, so that 5a = Sa&z- 

In all, we consider head-on collisions for seven differ- 
ent initial configurations of both equal (EQ) and unequal 
(NE) masses, and varied spins. The initial directly spec- 
ified parameters of all the simulations are listed in Table 
HI while derived quantities are given in Tablc[TTl all length 
scales are given in terms of a fiducial mass M [which co- 
incides with the Arnowitt-Deser-Misner (ADM) mass of 
the system for equal- mass runs] . For each run, the punc- 
tures were initially placed on the y axis with the center 
of mass at the origin. The proper separation I between 
the holes, along the initial slice, is measured between the 
closest parts of the apparent horizons along the y axis. 
We note that I is not necessarily the smallest physical 
distance between the holes (in this spatial slice), as spin 
effects can twist the space-like geodesies off the coordi- 
nate axis. 

We use the horizon mass m to define the black hole 
mass ratio q on the initial slice, 

q = m 1 /m 2 , (2) 



TABLE I: Directly specified parameters of the numerical sim- 
ulations, in terms of the fiducial mass M. S\ is the non-zero 
component of the Bowen-York angular momentum on each 
hole. NEa+_ and NEb+_ are, respectively, large-separation 

and equal-spin-parameter variants of NE-| , as explained in 

the text. 



Run 


mi >p /M m 2 , p /M y x /M 


?/ 2 /M 


Sf/M 2 Sf/M 2 


EQ 00 


0.5000 


0.5000 


4.0 


-4.0 


0.0 


0.0 


NEoo 


0.4909 


0.7478 


4.8348 


-3.2232 


0.0 


0.0 


EQ+o 


0.3444 


0.5000 


4.0 


-4.0 


0.2 


0.0 


EQ + _ 


0.3444 


0.3444 


4.0 


-4.0 


0.2 


-0.2 


NE+_ 


0.3436 


0.7140 


4.8 


-3.2 


0.2 


-0.2 


NEa+_ 


0.3436 


0.7140 


7.2 


-4.8 


0.2 


-0.2 


NEb+_ 


0.3436 


0.5496 


4.8 


-3.2 


0.2 


-0.4486 



where m\ is the horizon mass of the lighter hole. The 
symmetric mass ratio is 

v = mi m 2 /(mi + m 2 ) 2 ■ (3) 

The horizon mass is derived from the apparent hori- 
zon's irreducible mass using the Christodoulou formula 

s 2 

m 2 = m% r + — —, (4) 
4mf rr 



where m !rr = yj A ah /167T and A ah is the area of the ap- 
parent horizon, which we locate using an adapted version 
of Jonathan Thornburg's AHFinderDirect code [53| . 

The spin vector 5 applied to each hole in the Bowen- 
York data prescription {1} is in fact the ADM angular 
momentum of that hole - the total angular momentum 
of the spatial slice as measured by an ADM integral at 
infinity if no other sources were present. Thus these pa- 
rameters represent a global quantity. 

The more standard definition of astrophysical black- 
hole spin, however, is a local one. The Kerr solution is 
parametrized by the black hole horizon mass m, and a 
spin parameter a restricted to the range of values [0, to]. 
Then, a = |5|/to. We note, however, that it has proved 
impossible to bring Bowen-York spinning data to an a 
value of more than ~ 0.927m, significantly below the 
maximal Kerr value (54| . 

Even in the head-on case, with spin vectors orthogonal 
to the separation vector, there is much freedom in the 
parameters describing the initial data. In particular, we 
have chosen to scale all distances so that the horizon 
mass of the lighter hole, mi, is kept the same relative to 
the grid spacing (ensuring a common level of resolution of 
that hole's features). Within this restriction, we have also 
striven to maintain the same proper horizon-to-horizon 
separation I between the holes. 
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III. METHODOLOGY 



IV. CODE CALIBRATION 



These initial data sets were evolved using the moving 
puncture method as implemented in the Hahndol code 
[33 ]. We used standard Baumgartc-Shapiro-Shibata- 
Nakamura (BSSN) (5f| [56[ evolution equations, with the 
addition of dissipation terms as in [57] and constraint- 
damping terms as in [58j |. in order to ensure robust sta- 
bility. Our gauge conditions are 1+log lapse slicing and a 
hyperbolic Gamma-driver shift condition, as used in [44| . 
Here we take r\ = 2.0, and the initial lapse shape chosen is 
otinit = V'o" 4 , where ipo is the initial conformal factor; this 
gives a steeper fall-off near the punctures than the ip^ 2 
used in the evolutions of [32| . Time integration was car- 
ried out with a fourth-order Runge-Kutta algorithm, and 
spatial derivatives with fourth-order-accurate finite dif- 
ferencing stencils. We employed a second-order-accurate 
Sommerfeld condition on the outer boundary at 256A/; 
physical modes propagating in from the outer boundary 
will then take ~ 190M to reach the outermost extraction 
sphere at 60M, while gauge modes propagating at \/2 
times the speed of light will take ~ 140M to reach this 
same sphere. As will be seen, the bulk of the merger ra- 
diation and associated momentum flux will have passed 
through the extraction sphere by this time. We used 
adaptive mesh refinement (AMR) implemented via the 
software package PARAMESH [H, S E3, HI, with fifth- 
order-accurate interpolation between refinement regions 
(as preservation of fourth-order accuracy in the bulk de- 
mands better than fourth-order-accurate interpolation at 
refinement interfaces [46l.[63l.[6l]). 

The momentum kick of the merged black hole is calcu- 
lated as the aggregated time-integral of the momentum 
flux, or thrust, dP % jdt - a surface integral of the squared 
time-derivative of the radiative Weyl scalar ip^ times the 
unit radial vector [6511: 
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(5) 
(6) 



where to is the initial time in the simulation (all time 
integrals above should properly be from t' = —00; our 
implementation neglects thrust contributions from before 
to). In fact, for reasons discussed in the next section, 
we must delay the integration start time to to when the 
thrust passes through zero. As §5§ is quadratic in the 
waveform -04, the kick can be expected to be weak. To 
perform the angular integration in (J^Jl . we use the second- 
order Misner algorithm described in [66l . [67| ■ 



We have established the convergence of the Hahndol 
code in prior publications [32l . l3~ij . To these previously 
published results, we have added only one qualitatively 
new aspect: spin on the pre-merger holes. Thus we re- 
strict our discussion of convergence here to a study of 
the convergence properties of a sample spinning data set 
with equal masses and anti-aligncd spins, EQ^ . 

Our main evolutions were performed at a maximum 
resolution of M/32 in the vicinity of the pre-merger holes. 
We emphasize here that, as we have maintained a smaller 
horizon mass of m « 0.5 M in all runs, the grid spac- 
ing realized near the smaller hole will always be close 
to m/16. This yields corresponding grid spacings of 
h = M, h = 2M and h = 8M at the extraction surfaces 
-Roxt = 30M, i? cxt = 60M and the outer boundary, re- 
spectively. Relative to this standard resolution (M/32), 
we have added two higher resolutions (M/40 and M/48). 

Fig. [1] indicates the rate of convergence of the Hamil- 
tonian and momentum constraints in two regions of in- 
terest: the strong-field region surrounding the holes (left 
panels), and the refinement level containing our i? e xt = 
60M wave-extraction sphere (right panels) . In our simu- 
lations, the Hamiltonian constraint exhibits greater than 
third-order convergence, and the momentum constraint 
exhibits second-order convergence through most of the 
duration of the runs. This may be partly due to our 
second-order-accurate initial data solver; however errors 
of lower-than-expected order have also been attributed to 
"leakage" from the punctures [38, 68]. At late times the 
apparent convergence of the momentum constraint drops 
to a rate between first and second order; we believe this 
is because of a high-frequency gauge pulse propagating 
outward until it is insufficiently resolved at the grid res- 
olution achieved in the outer zones of our coarser run. 
However, this gauge pulse does not significantly affect 
the thrust computation, as evidenced by the invariancc 
with extraction radii demonstrated in Sec. IV1 

In Fig. [2[ we illustrate the convergence of the dom- 
inant (2,2) mode of the Weyl scalar Re(V>4), extracted 
at i? ex t = 60M. For comparisons, we have lined up the 
main amplitude peak in time across resolutions, and ro- 
tated the waveform phase to zero at this point (similar 
to the maquillage used in inspiral waveform comparisons 
[3ll]). Resolution-dependent errors are higher during the 
unphysical high-frequency Bowen-York pulse (see next 
section); by ~ 100 M, when the physical merger radia- 
tion reaches the detector, errors are consistent in overall 
magnitude with second-order convergence, and compare 
favorably with the equivalent plots in Fig. 3 of [30j . 



J 



4 



TABLE II: Derived initial parameters in terms of the fiducial mass M. mi is the horizon mass of the lighter hole, q is the 
mass ratio defined by horizon masses, oa = ISUI/wia is the approximate Kerr spin parameter of hole A. I is the proper 
horizon- horizon separation of the holes. Madm was calculated at a finite coordinate distance from the origin (6071/ for all but 
run EQ +0 ), using Eq. (12) of [5l[. Errors in horizon masses mi and nii are ~ 1.5% (see text), which propagate into derived 
quantities q, v, ai, 02, ai/mi, and ai/mi. ADM mass errors enter at the fourth significant digit. 



Run 


mi/M m 2 /M 


q 


V 


a\ 


2 


ai/mi 


(12/7712 


l/M 


Madm/M 


EQoo 


0.514 


0.514 


1.0 


0.250 


0.0 


0.0 


0.0 


0.0 


12.24 


1.0 


NEoo 


0.514 


0.771 


0.667 0.240 


0.0 


0.0 


0.0 


0.0 


12.24 


1.24 


EQ+o 


0.516 


0.514 


1.004 


0.250 


0.389 


0.0 


0.758 


0.0 


12.24 


1.0 


EQ+_ 


0.514 


0.514 


1.0 


0.250 


0.389 


0.389 


0.758 


0.758 


12.4 ± 0.2 


1.0 


NE+- 


0.516 


0.773 


0.668 


0.240 


0.388 


0.259 


0.752 


0.335 


12.6 ± 0.2 


1.24 


NEa+_ 


0.513 


0.764 


0.672 


0.240 


0.390 


0.262 


0.759 


0.342 


17.0 ± 0.2 


1.25 


NEb + _ 


0.516 


0.784 


0.658 


0.239 


0.388 


0.572 


0.752 


0.730 


13.0 ± 0.2 


1.26 



In Fig. [3l we illustrate the convergence of the mo- 
mentum thrust for our data set EQ_| , extracted at 

-Roxt = 60M (as for Fig. [2]); beginning at the integra- 
tion start time i = 90. 5M. To compare the different 
data sets, we re-zero and rescale the time axis to overlay 
the two largest physical peaks for all three resolutions. 
The differences have again been scaled for first- and 
second-order convergence. The thrust displays better 
than second-order convergence for most of the physical 
thrust, but has a significant excursion around t = 102717. 
Nevertheless, the relative errors in the thrust are ~ 1% 
(compare with upper panel). 



will mix with the physical radiation we are generally in- 
terested in. Nevertheless, since the pulse is released im- 
mediately, and is of short duration, it should be easily 
identifiable in the total waveform. 

To illustrate, we show in Fig. [4] the dominant (I = 
2, m = 2) spin-2-weighted spherical harmonic mode of 
the "outgoing" Weyl scalar ip4 (extracted at i? ox t = 
30M), for three different equal-mass configurations: 
zero-spin (EQ 00 ), single-spin (EQ +0 ) and double-spin 

(EQ_| ). Despite the differences in spins, these three 

waveforms agree well in the later "physical" part of the 
signal, with stronger differences occurring in the non- 
physical BY pulse. In the zero-spin case, the BY pulse is 
small, but it is significant for EQ +0 and EQ_| . 



A. Bowen-York Radiation Pulses 

Bowcn-York black hole initial data contains a non- 
negligible amount of "false radiation," even after the 
Hamiltonian constraint has been solved. This radiation 
is unphysical in the sense that it is a reflection of the 
unphysical approximations, such as conformal flatness, 
made in generating the initial data (c.f. [69|). 

To estimate the approximate duration of this unphys- 
ical radiation pulse, we assume that each hole relaxes to 
a Kerr hole (of approximately the same mass and spin) 
through the emission of quasinormal modes (QNMs). 
The most slowly damped of these will have an "e-folding" 
time T e « 12m for a hole of mass m [7(j • Then a pulse will 
have decayed by two orders of magnitude after ~ 55m - 
that is, ~ 30A7 for our equal-mass (EQ) data sets, and 
~ 45M for our unequal-mass (NE) data. 

This Bowen-York "pulse" is concentrated close to each 
hole's horizon, and much of it falls into its parent hole, re- 
sulting in an increase in apparent-horizon mass of ~ 1%; 
however, what radiation escapes to our detection spheres 



This same pulse will yield a non-zero contribution to 
the momentum kick. Again, this should be easy to iden- 
tify, as it will result in a plateau early in the run. Hav- 
ing identified a time to at which the BY pulse and its 
"pseudo-kick" are finished, we remove this effect by be- 
ginning the integration of |(5j) from that point. 



V. RESULTS 

We group our simulations into three categories. First 
we consider the effect of unequal masses alone, where 
the kick experienced should be along the line connecting 
the holes. We then add spins to holes of equal horizon 
mass, to explore the dependence of kicks on pure spin 
anisotropy; this configuration should yield kicks orthog- 
onal to the line connecting the holes. Finally we consider 
the combination of both spin and unequal masses, to ex- 
plore the interrelation of the two mechanisms generating 
kicks. 
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FIG. 1: Ll-norm convergence of Hamiltonian (top) and momentum (x component; bottom) constraints for the EQ H run, in 

two regions of the simulation grid: just outside the horizons (left), and in the primary radiation extraction region (right). The 
Hamiltonian constraint displays third-order convergence in both regions, while the momentum constraint convergence drops to 
1.5-order in the extraction region. 



A. Unequal masses without spin 

The head-on collision of two holes of unequal masses 
has been considered several times in the past, both with 
close-limit analysisfCLA) [7l[ and with fully numerical 
evolutions in 2D [72|]. The CLA results (performed with 
Misner two-sheeted data) provide a definite prediction of 
preferential radiation of linear momentum; the numeri- 
cal results seemed to confirm this, and indicated where 
the CLA fails. The numerical results were for proper 
separations far less than ours; extrapolating their results 
to the larger separations we treat here indicates that a 
longitudinal kick ~ lOkm/s would be produced. 

In Fig. [5j we plot the y-component of the thrust dP/dt 
from Eq. |(SJ), and its time integral from the NEoo run. 
Note that the x and z components of dP/dt are zero by 
symmetry: the kick is longitudinal, i.e. along the axis of 
collision (y axis). The final kick velocity, given in Table 
EH is ~ 2.7km/s. 



The final kick obtained from Fig. [5] is highly consis- 
tent across extraction radii, with a spread between the 
R cxt = 30M and R ext = 40M values of - 0.5%. We note 
however that the difference between i? cxt = 50M and 
-Rext = 60M is larger than that between R ext = 40M 
and i? cxt = 50M. This appears to be because the extrac- 
tion surface at i? ex t = 60M lies in a region of less grid 
refinement than the three closer surfaces. 



B. Equal masses with spin 

We now turn to simulations of spinning black holes. 

We ran two configurations (EQ +0 ,EQ H ), as indicated 

in Table |TJ As only one particle is spinning in the for- 
mer case, we should see only "spin-orbit" effects between 
that spin and the motion of the particle. In the latter 
case, with two spins, we might also see "spin-spin" ef- 
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FIG. 2: Convergence of Re(t/>4), the real part of the dominant 

(2, 2) mode of the waveform for the EQ H run. The upper 

panel shows the thrust for central resolution M/32, extracted 
at 7?ext = 60M; the lower panel shows the absolute differ- 
ences M/32 - M/40 (solid line) and M/40 - M/48, the latter 
scaled to compare with coarse-medium for 1st- and 2nd-order 
convergence. The inset focuses on the physical part of the 
waveform, which shows much smaller resolution errors. 
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FIG. 4: Dominant (I = 2, m = 2) mode of the waveform ^4 
for the three equal-mass configurations: EQ 00 (solid/black), 

EQ, (dashed/red), and EQ H (dot-dashed/green). These 

data were extracted at coordinate distance i? cxt = 30 M. 



PN order than for spin-orbit [731 ] 1 

In Fig. [S]we show the transverse kick AP X observed for 
both equal-mass data sets. We note that the magnitude 

of both the BY pulse and the final kick for the EQ H 

data are roughly double those for EQ +0 (see Table llllj) . 
This indicates that the total kick is close to a simple sum 
of the individual kicks. Thus any spin-spin contribution 
to the total kick is negligible in comparison to the spin- 
orbit term. 



C. Unequal masses with spin 

Having considered collisions between non-equal mass 
and non-spinning black holes and equal mass and spin- 
ning black holes, we combine both spin and mass ratio 
in a few cases, the NE_| , NEa_| , and NEb_| runs. 

We show in the upper panel of Fig. [7]the thrust for the 
NE H case. We note that unlike the simpler cases before, 



FIG. 3: Convergence of transverse momentum thrust dP x jdt 
for the EQ + _ run. The upper panel shows the thrust for 
central resolution M/32, extracted at i? ox t = 60M; the lower 
panel shows the absolute differences M/32 — M/40 (solid line) 
and M/40 — M/48, the latter scaled to compare with coarse- 
medium for 1st- and 2nd-order convergence. 



fects; however PN predictions for the latter arc of higher 



For the EQ_|_ run presented here, we have used slightly different 
choices in our numerical evolution: advection terms employed 
a third-order upwinding scheme, guardccll-filling was fourth- 
order-accurate, time stepping was performed with a second-order 
Crank-Nicholson method, no dissipation was used, and the ini- 
tial lapse function shape was cti n a = ^bL' as ' n ev °l u tions 
°f |32i ] . For this data alone, the physical outer boundary was at 
128M, but the outermost extraction here was at -R ex t = 30M. 
From direct comparison of the thrusts obtained from the two 

evolution schemes for the EQ_| data, we have determined that 

the only measurable difference is a gauge-driven IM delay in the 
arrival of wavefronts at the extraction sphere. 
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FIG. 5: Longitudinal thrust dP v /dt (top) and kick AP y (bot- 
tom) for NEoo, at extraction radii i? ox t of 30M, 40M, 50M 
and 60M, where the latter three have been time-shifted by 
10.6M, 21.0M and 31AM, respectively, to align the highest 
thrust peak with the 30M case (consistent with time-shift 
formula Eq. (14) in [13). 



the BY pulse and physical signal are not obviously sep- 
arated in the thrust. As a result, we may expect to have 
trouble determining where to begin the kick integrations. 

In principle, a larger initial separation would mean 
more time between the BY pulse - emitted at t = and 
lasting ~ 43M (see discussion in Section |IVA|) - and the 
bulk of the physical radiation, which only becomes large 
near merger. The main physical kick will certainly be 
different in magnitude between the two cases, and even 
the BY pulse magnitude may differ, depending on the 
dependence of the spectrum on the initial separation of 
the binary. Nevertheless, we can expect the BY pulse 
duration to be the same, as this is determined by the 
leading quasinormal frequencies of each hole, and these 
frequencies are fully determined by the mass and spin of 
each hole, which is a constant between the two cases. In 
Fig. [7J we show dP x /dt and AP X for NE H , with a coor- 
dinate separation of 8M, and NEa_| , with a coordinate 

separation of 12M. 

The extra initial separation (NEa^ ) leads to a very 

similar physical thrust, delayed ~ 22M relative to the 

original separation (NE_| ). This physical thrust begins 

at ~ 11AM for the detector at i? ex t = 60M. Integrating 
forward from any time up to ~ 20M before this will 

result in the same transverse kick for the NEa_| case. 

To ensure that as much of the physical kick as possible 
is obtained in the smaller-separation case, we integrate 
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-5e-05 ■ 
< 

-0.0001 




50 
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FIG. 6: Transverse thrust dP x /dt and kick AP X for the 
spinning, equal-mass cases EQ +0 (black/solid) and EQ + _ 
(red/dashed), at extraction radii 7? ex t of 30M, where the 
EQ + _ data have been time-shifted by 1M to compensate for 
gauge differences. 



the NE-| thrust from to 

i?ext = 60M. 



114M - 22M = 92M at 



In Figs. [8] and [9] we present the momentum kicks in 
the transverse (x) and longitudinal (y) directions for the 

NE-| data set, comparing with the kicks seen from the 

NE o an d EQ + _ cases before. 

Taking the transverse direction first, we see from Fig. 

|S] that the momentum kick from NE_j is significantly 

smaller than that of EQ_| . The spin angular momen- 
tum S present in each case is the same, indicating that 
the kick is not a function of S. In contrast, the momen- 
tum kick observed from NEb_| is much larger. 

The physical momentum kicks from these runs are pre- 
sented in Table IIIII We note that the ratio of trans- 
verse momentum kicks AP X between the cases NE_| and 

NEb-| is roughly 2/3, consistent with the PN-derived 

ratio (see Eq. (O below). 

In addition to finite differencing inaccuracy, for these 
cases we find that the high slope in the thrust for these 
cases yields a non-negligible error associated with the 
choice in integration starting point to. As before, we as- 
sess a timing-related error by integrating the thrust from 
to and from to ± 5M. The resulting central kicks and 
errors are given in Table IIIII 

With this in mind, we turn to the longitudinal direc- 
tion. Here we see that the initial BY burst is compara- 
ble in magnitude to the physical later part. Using for 
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FIG. 7: Transverse thrust dP x /dt (top) and kick AP X 
(bottom) for unequal-mass-with-spin data, at coordinate 
separations 8M (NE+_; black/solid) and 12M (NEa+_; 
red/dashed). These data were extracted at i? oxt = 60M. 



consistency the same integration starting time as in the 
transverse direction, we present the longitudinal kicks in 
Table Mil It is clear that the to uncertainty leads to 
large relative errors for the longitudinal kicks. To this 
uncertainty, all longitudinal kicks are consistent with the 
zero-spin case, NE o- 



D. Summary of Results 

We draw together in Table IIIII the kicks observed in 
each of our simulations. We can develop some expec- 
tations for the comparative results by referencing post- 
Newtonian theory estimates of the radiative linear mo- 
mentum loss due to unequal masses and spin-orbit cou- 
pling. For example, adapting Eqs. (3.31a-b) of [73| to 
radial infall along the y axis, we find: 



P 



N 



16 5m 2 ^n r iT\ 



105 rriT 



V r ) 



2niT 



yi 



Pso = - 



16 v 2 f 2 /tot\ 4 



15 r 



+ a 2 )e x , 



(7) 



(8) 



where mj- = mi+m 2 , 5m = mi—m 2 , and r is the spatial 
separation of the two particles, with r < 0. 

In studies of spin such as the one carried out here, 
whether physical effects scale with angular momentum 
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FIG. 8: Transverse thrust dP x /dt (top) and kick AP X (bot- 
tom) for NEoo (black), EQ H (red), and NE + _ (green), and 

NEb+- (blue). These data were extracted at Rext = 60 M. 



5 1 , the Kerr parameter a = S/m or the dimensionlcss 
number a/m is an open question. The post-Newtonian 
result of Eq. predicts a thrust (and hence kick) that 
scales with a for each hole. We have tried to address this 
uncertainty in our choice of data sets. 

It is informative to compare the transverse kick re- 
sults from each of our spinning data sets with each other, 
in light of post-Newtonian predictions. Looking at Eq. 
©, we assume that the right-hand side is more or less 
insensitive to the details of the rate of infall, the im- 
portant length scale then being the total mass of 
the binary. With this in mind, the remaining freedom 
lies in the spin scaling (ai + a-i)jmT. Relative to the 

EQ_| case, this scaling factor is {1/2,2/3,2/3,1} for 

cases {EQ +0 ,NE + _,NEa + _,NEb+_}. In Fig. [TO] we 
combine the post-BY-pulse part of the transverse thrust 
dP x /dt of each case, with times rescaled according to the 
ADM mass from Table[H and translated so that the peaks 
coincide. We have also rescaled the thrust amplitude by 
the factors appropriate to the PN above Surprisingly, af- 
ter this rescaling, the five thrusts fit very well, with a 
deviation of overall amplitude of ~ 20%. This appears to 
indicate that PN predictions of recoil kicks have validity 
further into the merger regime than would be expected. 



We have looked at the puncture tracks of the holes, 
which give a coordinate speed estimate for the final hole. 
These agree with the radiation-generated kick in order- 
of-magnitude and relative scaling (transverse vs longitu- 
dinal velocity components), but not in detail; in partic- 
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ular, we do not see a velocity plateau; given the small 




recoil velocities (compared to inspiral) in question here, 
the puncture-derived position may just not be accurate 
enough. It should be borne in mind also that the punc- 
ture track is highly coordinate-dependent, and has no 
guaranteed relationship to physical quantities such as the 
kick. 



t/M 



FIG. 9: Longitudinal thrust dP y /dt (top) and kick AP y (bot- 
tom) for NEoo (black), EQ H (red), and NE+- (green), and 

NEb+_ (blue). These data were extracted at J? ox t = 60 M. 
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FIG. 10: Transverse thrust dP x /dt for all spinning data sets: 
EQ +0 , EQ + _, NE+_, NEa + _, and NEb+_. Time has been 
rescaled by the ADM mass for each data set, and translated 

relative to EQ_j to line up peaks. Amplitudes have been 

scaled relative to EQ_| according to PN predictions. These 

data were extracted at i? oxt = 30 M . 



We also present in Table flVl estimates for the final Kerr 
parameters Mf and aj of the post-merger hole for each 
simulation. We calculate these estimates in two differ- 
ent ways: First, we compare the initial ADM energy and 
angular momentum of the system with the radiated en- 
ergy and angular momentum. In this we arc limited by 
the accuracy of the ADM energy, reliable to three signif- 
icant digits. Second, we use the (2, 2) mode of the radi- 
ation after the final amplitude peak to extract a quasi- 
normal-mode frequency and damping time, and hence 
to determine the underlying Kerr parameters, assuming 
the perturbative limit. This latter method is sensitive 
to the determination of an appropriate starting time for 
the QNM fit, as well as the neglect of QNM overtones. 
The determination a/M — 0.0252 for the NEoo data set, 
which can contain no angular momentum, indicates the 
level of error of this method. 



J 



VI. DISCUSSION 



In this paper, we have addressed the roles of unequal 
masses and spins in producing recoil kicks in head-on col- 



lisions. We have chosen sample spins and mass ratios to 
explore how important spin effects are for kicks, and how 
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TABLE III: Final integrated momentum kicks and corresponding kick velocities. In each case, we have removed "BY pulse" 
effects through deferring integration until after the passage of the pulse. Quoted errors at each extraction radius are obtained 
by varying the integration starting time by ±5M. 



Run 


to 
(M) 


AP X AP y 
(10" 5 ) (10" 5 ) 


V X Vy V 

(km/s) (km/s) (km/s) 


NE 00 30M 
40M 
60M 


51.0 
61.7 
82.9 


0.0 1.12 ±0.02 
0.0 1.14 ±0.04 
0.0 1.16 ±0.04 


0.0 2.71 ±0.05 2.71 ±0.05 
0.0 2.76 ±0.10 2.76 ±0.10 
0.0 2.81 ±0.10 2.81 ±0.10 


EQ +0 20M 
30M 


46.1 
56.6 


-5.13 ±0.04 0.0 
-5.51 ±0.01 0.0 


-15.38 ±0.12 0.0 15.38 ±0.12 
-16.52 ±0.03 0.0 16.52 ±0.03 


EQ + _ 30M 
40M 
60M 


59.6 
70.0 
90.5 


-10.70 ± 0.08 0.0 
-10.53 ± 0.09 0.0 
-10.81 ±0.15 0.0 


-32.08 ±0.24 0.0 32.08 ± 0.24 
-31.57 ±0.27 0.0 31.57 ±0.27 
-32.41 ±0.45 0.0 32.41 ±0.45 


NE+_ 40M 
60A/ 


72.0 
92.0 


-8.47 ±0.31 1.36 ±1.10 
-8.50 ±0.27 1.52 ±1.36 


-20.41 ±0.75 3.28 ±2.65 20.67 ± 0.85 
-20.48 ±0.67 3.66 ±3.28 20.80 ± 0.88 


NEa+_ 40M 
60M 


94.0 
114.0 


-7.89 ±0.26 1.04 ±0.08 
-7.89 ±0.21 1.00 ±0.07 


-18.92 ±0.63 2.49 ±0.20 19.08 ± 0.63 
-18.92 ±0.51 2.40 ±0.17 19.07 ± 0.51 


NEb+_ 40M 
60M 


72.0 
92.0 


-12.33 ± 0.56 2.03 ±0.64 
-12.44 ± 0.39 2.00 ±0.55 


-29.41 ± 1.34 4.84 ± 1.53 29.81 ± 1.35 
-29.67 ±0.93 4.77 ± 1.31 30.05 ± 0.94 



TABLE IV: Basic parameters Mj and o/ of the post-merger hole, as calculated from radiated energy and angular momentum 
("rad"), and from measured quasi-normal mode frequencies ("QNM"). All quantities based on extraction at 7? e xt = 60Af, 
except for _R ex t = 30M for run EQ +0 . 



Run 


Madm/M 


Erad/M 




M f , rad /M Jadm/M 2 


Jrad/M 2 




(J/M 2 ) f , rad 


M f}QNM /M (a/M)f, QNM 


NEoo 


1.24 


-6.14 x 10" 


-4 


1.24 


0.0 


0.0 




0.0 


1.2536 


0.0252 


EQ+o 


1.00 


-7.90 x 10" 


-4 


1.00 


0.2 


-1.101 x 10 


-3 


0.1992 


1.0102 


0.3507 


EQ+_ 


1.00 


-9.02 x 10" 


-4 


1.00 


0.0 


-0.000 x io- 


-4 


0.0 


1.0113 


0.0198 


NE+- 


1.24 


-8.19 x 10" 


-4 


1.24 


0.0 


-3.508 x 10' 


-4 


2.33 x 10~ 4 


1.2585 


0.0261 


NEa+_ 


1.25 


-8.16 x 10" 


-4 


1.25 


0.0 


-1.582 x 10 _ 


-4 


1.01 x 10" 4 


1.2676 


0.0300 


NEb+_ 


1.26 


-1.10 x 10" 


-3 


1.26 


0.2486 


-5.000 x 10" 


-4 


0.154 


1.2582 


0.2175 



the effects of spin and mass ratio are related in gener- 
ating kicks. We have carried out kick extractions from 
data sets with mass ratios of 1 and 2/3, and with spins 
on one or both holes. 

We have observed that head-on collisions of spinning 
black holes can produce significant kicks. For the anti- 
aligned spin cases we have studied, the spin-kicks are 
transverse to the direction of initial separation. The mag- 
nitude of the kicks, and even the thrust curves (Fig. [TO)) , 
scale with the sum of the individual black hole spin pa- 
rameters, ai ± 02. These kicks easily exceed the longi- 
tudinal kick produced in the case of unequal masses. To 
the accuracy available in our simulations, we find that 
the kicks due to the mass ratio and spins are indeed or- 
thogonal and independent. Both the spin scaling and 
the orthogonality of spin-kicks and mass-ratio-kicks are 
consistent with leading-order PN predictions for these ef- 
fects, as given by Eqs. I|7I8[) . 

It is remarkable that the PN predictions seem to de- 
scribe our results. The kick-producing radiation derived 



in our simulations is generated in the systems' non-linear 
mergers and ringdowns, where the assumptions behind 
the PN approximations clearly do not apply. In gen- 
eral terms, the PN analysis (see [73|]) provides that the 
unequal-mass-kick is produced by a coupling of the mass- 
quadrupole and mass-octupole moments, while the spin- 
kick is produced by a coupling of the mass-quadrupolc 
and current-quadrupole moments. At least heuristically, 
we can consider these moments even in the nonlinear 
problem. In our problem, the presence of spin does 
not seem to have a large effect on the mass moments, 
as is suggested by the spin-independence of the mass- 
quadrupole-dominated (I = 2,m = 2) waveforms, seen 
in FigfJ] The leading effect of putting spin on the black 
holes is to scale up the current-quadrupole moment. For 
general spins, this term scales with {S\/m\ — S^/""^)- 

In the more important case of inspiraling black holes, 
our results suggest that the current-quadrupole effects 
should also provide strong kicks, scaling in a similar way 
with spin. Because of the scaling with S/m we would 
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generally expect the strongest spin-kicks for nearly-equal- 
mass mergers. Our results support the use of PN expres- 
sions (as in [l^,[4l|,|42|) to predict the scaling of spin- kicks 
for black hole configurations that have not yet been stud- 
ied. 

Since initial submission of this paper, several groups 
have begun mode-based analysis of recoil kicks from bi- 
nary inspiral [zH, [75| . Additionally, further investigation 
of the maximum size and generi city of spin supcrkicks 
has been undertaken [76|, O, [ll, . These substantial 
predictions have opened the door to exciting new possi- 
bilities in the astronomy of supermassive black holes and 
galactic nuclei. 
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